Seismic Rayleigh waves on an 
exponentially graded, orthotropic half-space 



Michel Destrade 
2006 

Abstract 

Efforts at modelling the propagation of seismic waves in half-spaces 
with continuously varying properties have been mostly focused on 
shear-horizontal waves. Here a sagittaly polarized (Rayleigh type) 
wave travels along a symmetry axis (and is attenuated along another) 
of an orthotropic material with stiffnesses and mass density varying 
in the same exponential manner with depth. Contrary to what could 
be expected at first sight, the analysis is very similar to that of the 
homogeneous half-space, with the main and capital difference that the 
Rayleigh wave is now dispersive. The results are illustrated numer- 
ically for (i) an orthotropic half-space typical of horizontally layered 
and vertically fractured shales and (ii) for an isotropic half-space made 
of silica. In both examples, the wave travels at a slower speed and pen- 
etrates deeper than in the homogeneous case; in the second example, 
the inhomogeneity can force the wave amplitude to oscillate as well as 
decay with depth, in marked contrast with the homogeneous isotropic 
general case. 

1 Introduction 

Love (1911) showed that a inhomogeneous half-space, consisting of an elas- 
tic layer covering a semi-infinite body made of a different elastic material, 
can sustain the propagation of a linearly polarized (shear horizontal) sur- 
face wave. The Love wave is faster than the elliptically polarized (vertical) 
Rayleigh (1885) wave and it has been observed countless times during earth- 
quakes or underground explosions. Another recorded phenomenon is that 
Rayleigh waves are dispersive, a characteristic which is incompatible with the 
context of a homogeneous half-space given by Rayleigh (1885): Love showed 



that his layer /substrate configuration could also support a two-partial, ver- 
tically polarized, surface wave. Because this configuration introduces a new 
characteristic length, the layer thickness h (say), a dispersion parameter is 
now kh where A; is the wave number, and that surface wave is dispersive. 

Subsequent analyses introduced more and more layers to refine the model, 
until is was considered practical to view the inhomogeneity of the half-space 
as a continuous variation of the material properties (Ewing et al. 1957). 
Chief among these continuous variations is the one for which the elastic 
stiffnesses and the mass density vary exponentially with depth, all in the 
same manner, proportional to a common factor exp(— 2aa;2) say, where a is 
the inverse of a inhomogeneity characteristic length, and X2 is the coordinate 
along the normal to the free surface, so that here a dispersion parameter 
is now a/k for instance. Hence Wilson (1942), Deresiewicz (1962), Dutta 
(1963), Bhattacharya (1970), and many others studied the propagation of 
surface waves in such inhomogeneous media; they were however interested in 
shear-horizontal waves (Love-type). The literature on Rayleigh-type surface 
waves in that type of media is quite scarce, probably because the difficulty 
exposed below is encountered quite early in the analysis. 

In an anisotropic elastic body with continuously variable properties, the 
general equations of motion read 

CijklUi^kl + CijkljUi^k = PUi^tt, (1) 

where u is the mechanical displacement, and Cijki and p are the elastic 
stiffnesses and the mass density, respectively. Now consider the propagation 
of an inhomogeneous plane wave with speed v and wave number k in the 
Xi-direction, and with attenuation in the a:2-direction, 

u = f/V'=(^i+5^2-^*), (2) 

in a half-space X2 ^ made of an orthotropic0 material with an exponential 
depth profile, 

{cii(x2), 022(2:2), C12 (0:2), ceQ{x2),pix2)} = e'^^'^^jcn, C22, 0^2, Cgg, P°}- (3) 

Here the Xi, X2, x^ directions are aligned with the axes of symmetry, a is a 
real number, and the c°j and p° are constants; also, U° is a constant vector 
and q a complex number so that the attenuation factor is k'^{q). Then the 
equations of motion ([T]) yield 

+ c?i - pv^ + 2i{a/k)qcl, q{cl^ + c°J + 2i{a/k)cl^ ] U° = 
g(c?2 + c^e) + 2iia/k)cl2 c°22q' + 4, - pv' + 2i{a / k)qc°22\ 

(4) 



""^An anisotropic material belongs to the orthotropic symmetry class when it possesses 
three mutually orthogonal planes of mirror symmetry. 
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At a = 0, the material is homogeneous, and the associated determinantal 
equation - the propagation condition - is a real quadratic in q"^ which can be 
solved exactly (Sveklo, 1948). 

At a 7^ 0, the propagation condition is seemingly a quartic in q with com- 
plex coefficients, whose analytical resolution might appear to be a daunting 
task and to preclude further progress toward the completion of a boundary 
value problem (note that it remains a quartic even when the material is iso- 
tropic.) Hence, Das et al. (1992) and Pal & Acharya (1998) stopped their 
analytical study of that problem at that very point. In fact the transfor- 
mation of the quartic to its canonical form reveals that it is a quadratic in 
q + i(a/fc), with real coefficients. That this is so has rarely been identified: 
Biot (1965), in the context of incremental static deformations, seems to be 
the only one who has recognized this simplification. The present paper shows 
that the Stroh (1962) formulation of this problem, combined with a change 
of unknown functions, leads naturally to the biquadratic in question. Then 
the propagation condition can be solved exactly, and the general solution 
of form ([2]) to the equations of motion follows. In particular, the resolution 
of the dispersive Rayleigh wave boundary value problem poses no particular 
difficulty after all. Section 2 exposes this analysis, and Section 3 applies it 
to two types of exponentially graded half-spaces: one which would be made 
of orthotropic shales if a ^ and another which would be made of silica 
(isotropic). There, it is seen for both examples that the infiuence of the in- 
homogeneity is more marked upon the wave speed (rapidly decreasing with 
a/k) than upon the attenuation factors (slowing increasing with a/k). It is 
also found that the attenuation factors for the displacement amplitudes are 
distinct from those for the traction amplitudes, and that the amplitudes can 
decay in an oscillating manner for the isotropic silica. These two features 
are unusual and are clearly due to the inhomogeneity. 

The overall aim of the paper is to show that simple, analytical, exact re- 
sults can be obtained for seismic Rayleigh wave propagation in an anisotropic, 
inhomogeneous Earth. Of course it is unlikely any "real" inhomogeneity can 
be such that the stiffnesses and the mass density all vary in the same man- 
ner as in dn]), because it then leads to bulk wave speeds (proportional to 
the square root of stiffnesses divided by the density) which are constant 
with depth. The analysis of more realistic models must turn to numerical 
simulations such as those based on the finite difference technique or on the 
pseudospectral technique or on techniques with Fourier or other function ex- 
pansions (e.g. Tessmer 1995). These methods however encounter difficulties 
for the implementation of accurate boundary conditions and of strong het- 
erogeneity. The spectral element method seem to alleviate those difficulties 
but, as stressed by Komatitsch & al. (2000), it must be validated against 
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analytical solutions. Such a solution validation procedure is indeed a cru- 
cial necessity of numerical simulations in geophysics, where different software 
packages can give widely different predictions (Hatton 1997). 



2 The dispersion equation 

Consider the propagation of a Rayleigh wave, traveling with speed v and wave 
number k in the xi-direction, in an inhomogeneous half-space X2 ^ made 
of the orthotropic material presented in the Introduction. The associated 
mechanical quantities are the displacement components Uj and the traction 
components (Jj2 {j = 1,2). They are now taken in the form 



(5) 



where the Uj, tj2 (j = 1,2) are yet unknown functions of X2 alone, to be 
determined from the equations of motion and from the boundary conditions. 

The equations of motion: aijj = pui^tt, can be written as the second-order 
differential system ([T]), or as the following first-order differential system. 



U'' 
t' 



kNi e^"^2jV2 
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(6) 



Here Ni, N2, K are the usual constant matrices of Stroh (1962), given by 



-1 



^12 



-22 



No 



- 1 







1 


r° 
^66 







•-22-1 



K 



X 



(7) 



„o2 
•-12 



where c° : = — and X 

as 



^22 

^(^2) 
becomes 



p°v^ . With the new vector function ^, defined 



[e-"^^C/(x2),e"^^t(a;2)]*, 



the system 

^' = \kN^ where 



X :-- 



Ni + i{a/k)I {l/k)N2 
kK Nl - i{a/k)I 



(8) 



(9) 



Hence the apparently anodyne change of unknown functions ([H]) transforms 
the differential system with variable coefficients (E]) into one with constant 
coefficients. 
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Now solve the differential system 
cent form, 



with a solution in exponential evanes- 



^{X2) = e^'^^X, ^{p)>\a\/k, 



(10) 



where C is a constant vector, p is a scalar, and the inequahty ensures that 



u{oo) = 0, t{oo) = 0, |(oo) = 0, 



'11^ 



because by (IHl) and ( fTOl) i . u{x2) behaves as: expfc(ip + a/k)x2 and t{x2) 
behaves as: exp k{ip — a/ k)x2- Note in passing that, in sharp contrast to the 
homogeneous case, the displacement field and the traction field have different 
attenuation factors: for u it is k[Q{p) — a/k]; for t it is k[Q{p) + a/k]. 

Then ^ and p are solutions to the eigenvalue problem: A^<^ = p(^. The 
associated determinantal equation is the propagation condition, here a bi- 
quadratic (and not a quartic as Eq.(jll) suggested). 



P 



Sp^ + P = 0, 



(12) 



where 



P 



+ cl,)X]/{cl,cl,)-2{a/k)\ 

— {a/k) [c°2 — 2C^2'^66 ~ '^ll'^22 + ('^22 + '^Gg)"^]/ ('^22'^66) 

+ {a/kf. 



(13) 



Let pi and p2 be the two roots of (|T2|) satisfying inequality ( ITOl) . That pair 
may be in one of the two forms: pi = ibi, p2 = ib2, or pi = —a+ib, p2 = a+ib, 
where b, bi, 62 are positive. In both cases, piP2 is a real negative number and 
Pi +P2 is a purely imaginary number with positive imaginary part. It follows 
in turn that 



P1P2 



.2^2 



pfP2 = -VP, Pi+P2 = W-iPl + Pi? = iy 2Vp - S. (14) 



The associated eigenvectors C^, are determined from: N(^^ = PjC'' , as 



p| + 2i{a/k)pj - Co 
-[P^j + K<^/k)pj + fipj + i{a/k)fo] 
-k[9iPj + K'^/k)go] 
-k[Xp] + /lo] 



(15) 
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where the non-dimensional quantities Cq, /i, /o appearing in the displacement 
components are given by 

e, = {a/kf + cl,{cl,-X)/{cl,cl,), 



f, = {a/ky + ic^-X)/cl,~cycl 
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f, = {a/kf + {c^-X)/cl, + clJcl,, (16) 

and the quantities gi, Qq, ho (dimensions of a stiffness) appearing in the 
traction components are given by 

^?l = C°-(l + C°2/cyX, 
5(0 = C° - (1 - C12/C°22)X, 

ho = {a/kyx - (c° - X)(c°6 - X)/cl,. (17) 
Now construct the general solution to the equations of motion ([9]) as 

€(^2) = 7ie*'^^"^C' + 72e*'^^"^C', (18) 

where the constants 71, 72 are such that the surface X2 = is free of trac- 
tions: t(0) = or equivalently: ^(0) = [C/(0),0]*. This condition leads to 
a homogeneous linear system of two equations for the two constants, whose 
determinant must be zero. After factorization and use of (JT4l) . the dispersion 
equation follows as 



giixVP + ho) + {a/k)goX\/2VP -S = 0. (19) 

This equation is fully explicit (X is the sole unknown) because P and 5* 
are given in ffT^ and gi, go, ho are given in ffT7|) . and it is clearly dispersive 
due to the multiple appearance of the dispersion parameter a/k. At a = 
(homogeneous substrate), it simplifies to 




X) (cl, -X)_ _ (c°-X)(cg6-X) ^ ^20) 



22 "-ee "^66 



the classic (non-dispersive) secular equation for Rayleigh waves in orthotropic 
solids. 



3 Examples: exponentially graded shales and 
silica 

As two examples of application, consider in turn that the half-space is made of 
a material with exponentially variable properties which is (i) with orthotropic 
symmetry and (ii) isotropic. 
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In Example (i) the starting point is a model proposed by Schoenberg and 
Helbig (1997), accounting for the vertical fine stratification and the vertical 
fractures found in many shales. In their numerical simulations, they used 
the following orthotropic elastic stiffness matrix. 
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(21) 



Note that here the matrix is density-normalized so that its components have 
the dimensions of squared speeds, expressed in (km/s)^ (Schoenberg 1994). 
Schoenberg and Helbig remark that "the rock mass behaves as if it contains 
systems of parallel fractures increasing the compliance in some directions"; 
integrating this information, a is assumed positive here. Also, (12T!) is as- 
sumed to be the elastic stiffness matrix on the free surface X2 = 0. 

In Example (ii), the half-space is assumed to be made of an exponentially 
graded material such that at the boundary, c^i = 7.85, c^2 = 1-61 (10^° N/m^) 
and p° = 2203 kg/m^ as in silica (Royer & Dieulesaint 2000). Here too, a is 
taken positive. 

If the half-spaces were homogeneous, then the Rayleigh wave would travel 
with speed v° = Xj p° where X is given by fpU|) . that is v° = 1.412 km/s 
for shales and v° = 3409 m/s for silica. For any given dispersion parameter 
a/k, the dispersion equation ( |T9l) in the inhomogeneous half-spaces gives a 
unique root X. In both examples, it has then been checked that for that 
X, the propagation condition (|T^ gives two roots such that the inequality 
(|TU|) 9 is always satisfied. Thus the surface wave exists for arbitrary value of 
a/k, and it travels with speed v = ^ Xj p°. Although this state of affair 
is acceptable mathematically, it seems reasonable to limit the range of a/k 
to values where the wave amplitude decays faster than the inhomogeneity. 
Because the amplitudes of the tractions t decay as exp —k[Q{p) +a/k], they 
always decrease faster than exp —2ax2 by ffTUl) 9: on the other hand, the 
amplitudes of the displacements u decay as exp —k[Q{p) — a/k]: thus they 
decrease faster than the inhomogeneity as long as Q{p) > 3a/ k. In Example 
(i), it turns out that this latter inequality is verified for a/k < 0.107, and in 
Example (ii), for a/k < 0.274. 

Fig. 1 shows the variation of the wave speed (decreasing) and of Q{pi), 
3(^2) (increasing) in Example (i) over the range O^a/fc^O.l. It has also 
been checked there that the attenuation factors for both the displacements 
amplitudes {k[Q{p) — a/k]) and the tractions amplitudes {k[Q{p) — a/k]) 
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increase also. In conclusion, the surface wave travels at a slower speed in the 
inhomogeneous shales than in the homogeneous shales, and it is less localized. 

Fig. 2 shows the variation of the wave speed (decreasing) and of Q{pi), 
53(^2) (increasing) in Example (ii) over the range ^ a/A; ^ 0.25. It has 
been checked again that the attenuation factors for both the displacements 
amplitudes {k['^{p) — a/k]) and the tractions amplitudes {k[Q{p) — a/k]) 
increase also. Here again, the surface wave travels at a noticeably slower 
speed in the inhomogeneous case than in the homogeneous case, and it is 
slightly less localized. A most interesting phenomenon occurs aX a/k ^ 0.211 
where the nature of the roots changes from the form: pi = ibi, p2 = ib2, to the 
form: pi = —a+ib, p2 = a+ib, so that the amplitudes switch from decaying in 
a real exponential manner to decaying in an exponential oscillating manner. 
This latter situation never arises in a homogeneous isotropic half-space. 
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Figure 1: Exponentially graded orthotropic shales: variations with the dis- 
persion parameter a/k oi (a) the surface wave speed and (b) the imaginary 
parts of the quantities pi and p2 appearing in (|T8l) (the dashed line is the plot 
of 3a/ k, above which Q{pi), '^{p2) must be for the wave to decrease faster 
than the inhomogeneity) . 




Figure 2: Exponentially graded silica: variations with the dispersion param- 
eter a/k of (a) the surface wave speed and (b) the imaginary parts of the 
quantities pi and p2 appearing in (|T^ (the dashed line is the plot of 2>a/k). 



10 



